% Appendices
% This file is included by bayesian_vs_3dvar_cn_main.tex

\newpage
\appendix

\section{FLEXINVERT 系统矩阵构造详解}
\label{app:flexinvert-matrices}

本附录详细介绍 FLEXINVERT 系统中三个关键矩阵的构造方法，包括数学公式、算法实现和计算复杂度分析。

\subsection{先验误差协方差矩阵 $\mathbf{B}$ 的构造}

\subsubsection{数学结构}

先验协方差矩阵 $\mathbf{B} \in \mathbb{R}^{N \times N}$ 表示通量先验的不确定性与空间-时间相关性，其数学结构为
\begin{equation}
B_{ij} = \sigma_i \sigma_j \rho_{ij},
\end{equation}
其中 $\sigma_i, \sigma_j$ 为通量误差标准差，$\rho_{ij}$ 为相关系数。

\paragraph{1. 通量误差标准差}

采用相对误差模型：
\begin{equation}
\sigma_i = \max(\alpha \cdot |x_{b,i}|, \sigma_{\min}),
\end{equation}
其中：
\begin{itemize}
\item $\alpha$：相对误差（如 $\alpha = 0.5$ 表示 50\% 不确定性）
\item $x_{b,i}$：先验通量值
\item $\sigma_{\min}$：最小绝对误差下限，防止零先验导致零方差
\end{itemize}

{\color{red}\textbf{[相对误差模型的合理性]：通量误差通常与通量大小成正比——强源区（如热带雨林）的绝对误差大于弱源区（如沙漠），但相对误差可能相近。$\sigma_{\min}$ 下限确保即使先验为零的网格（如极地冰川）也有非零误差，允许反演探索这些区域。}}

\paragraph{2. 空间相关函数}

采用指数衰减相关模型，考虑海陆差异：
\begin{equation}
\rho_{\text{spatial}}(i,j) = \exp\left(-\frac{d_{ij}}{\sigma_{\text{eff}}}\right),
\end{equation}
其中 $d_{ij}$ 为网格 $i, j$ 间的球面距离（km），有效相关长度为
\begin{equation}
\sigma_{\text{eff}} = \begin{cases}
\sigma_{\text{land}} & \text{同为陆地}（\sigma_{\text{land}} \approx 500\text{--}2000\,\text{km}） \\
\sigma_{\text{ocean}} & \text{同为海洋}（\sigma_{\text{ocean}} \approx 1000\text{--}5000\,\text{km}） \\
\sqrt{\sigma_{\text{land}} \sigma_{\text{ocean}}} & \text{陆海交界}（\text{几何平均}）
\end{cases}
\end{equation}

球面距离通过 Haversine 公式计算：
\begin{equation}
d = 2R \arcsin\left(\sqrt{\sin^2\left(\frac{\Delta\phi}{2}\right) + \cos\phi_1 \cos\phi_2 \sin^2\left(\frac{\Delta\lambda}{2}\right)}\right),
\end{equation}
其中 $R = 6371\,\text{km}$ 为地球半径，$\phi, \lambda$ 为纬度、经度。

{\color{red}\textbf{[空间相关的物理意义]：陆地通量误差的相关长度通常小于海洋，因为陆地生态系统异质性更高（森林、农田、城市的拼贴）。海洋通量由大尺度海洋环流和温度场控制，相关尺度更大。陆海交界使用几何平均避免过强相关，反映不同表面类型的相对独立性。}}

\paragraph{3. 时间相关函数}

采用指数时间衰减：
\begin{equation}
\rho_{\text{temporal}}(t_1, t_2) = \exp\left(-\frac{|\Delta t|}{\tau}\right),
\end{equation}
其中 $\Delta t = |t_1 - t_2|$ 为时间间隔（天），$\tau$ 为时间相关尺度。典型值：
\begin{itemize}
\item 陆地生态系统通量：$\tau \approx 7$--30 天（受生长季节、物候变化控制）
\item 海洋通量：$\tau \approx 30$--90 天（受海温缓慢变化控制）
\item 人为排放：$\tau \approx 1$--7 天（受经济活动周期控制）
\end{itemize}

{\color{red}\textbf{[时间相关尺度的选择]：$\tau$ 反映通量变化的``记忆''时间。植被光合作用受多日累积降水和温度影响，$\tau \sim 7$--14 天合理。但对于化石燃料燃烧，工作日与周末模式差异显著，$\tau < 7$ 天更适合。实践中可通过变分法优化 $\tau$，或从高频塔通量数据估计。}}

\paragraph{4. 完整协方差矩阵构造}

对于时空状态向量 $\mathbf{x} \in \mathbb{R}^{N_{\text{box}} \times N_t}$，完整协方差矩阵为：
\begin{equation}
B_{(t_1,i),(t_2,j)} = \sigma_{t_1,i} \sigma_{t_2,j} \cdot \rho_{\text{spatial}}(i,j) \cdot \rho_{\text{temporal}}(t_1, t_2).
\end{equation}

这是 Kronecker 积结构的推广：
\begin{equation}
\mathbf{B} = \mathbf{C}_{\text{temporal}} \otimes \mathbf{C}_{\text{spatial}},
\end{equation}
但加入了逐网格的误差标准差对角矩阵。

\paragraph{5. 本征值分解与存储优化}

完整协方差矩阵规模巨大（$N \sim 10^4$--$10^6$），FLEXINVERT 采用本征值分解压缩存储：
\begin{equation}
\mathbf{B} = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^T \approx \sum_{j=1}^{K} \lambda_j \mathbf{u}_j \mathbf{u}_j^T,
\end{equation}
仅保留满足 $\lambda_j > \epsilon \cdot \max(\lambda)$ 的 $K$ 个主模态（$\epsilon \sim 10^{-3}$--$10^{-6}$）。

矩阵-向量运算通过本征空间实现：
\begin{align}
\mathbf{B}^{1/2}\mathbf{v} &= \sum_{j=1}^{K} \sqrt{\lambda_j} (\mathbf{u}_j^T \mathbf{v}) \mathbf{u}_j, \\
\mathbf{B}^{-1/2}\mathbf{v} &= \sum_{j=1}^{K} \frac{1}{\sqrt{\lambda_j}} (\mathbf{u}_j^T \mathbf{v}) \mathbf{u}_j.
\end{align}

\paragraph{6. 本征值截断准则}

确定保留模态数 $K$ 的方法：
\begin{enumerate}
\item \textbf{方差阈值法}：保留累积方差占比 $\geq \theta$ 的模态（$\theta \approx 0.99$）
\begin{equation}
\sum_{j=1}^{K} \lambda_j \geq \theta \sum_{j=1}^{N} \lambda_j.
\end{equation}
\item \textbf{相对本征值法}：保留 $\lambda_j > \epsilon \max(\lambda)$ 的模态（$\epsilon \sim 10^{-6}$）
\item \textbf{固定秩法}：直接指定 $K$（如 $K = \min(M, \lfloor 0.1N \rfloor)$）
\end{enumerate}

{\color{red}\textbf{[本征值截断的优势]：(1) 存储从 $\mathcal{O}(N^2)$ 降至 $\mathcal{O}(KN)$；(2) 消除近零本征值，改善条件数从 $\kappa(\mathbf{B}) \sim 10^{6}$ 降至 $\kappa(\mathbf{B}_K) \sim 10^{3}$；(3) 矩阵-向量乘积复杂度从 $\mathcal{O}(N^2)$ 降至 $\mathcal{O}(KN)$。对于 $N=10^5, K=10^3$，存储减少 99\%，计算加速 100 倍。}}

\subsection{观测误差协方差矩阵 $\mathbf{R}$ 的构造}

FLEXINVERT 假设观测误差无相关，采用对角矩阵：
\begin{equation}
\mathbf{R} = \text{diag}(\sigma_1^2, \sigma_2^2, \ldots, \sigma_M^2).
\end{equation}

\subsubsection{误差分量}

观测误差包含多个分量：
\begin{equation}
\sigma_i^2 = \sigma_{\text{instrument}}^2 + \sigma_{\text{representation}}^2 + \sigma_{\text{transport}}^2,
\end{equation}
其中：
\begin{itemize}
\item $\sigma_{\text{instrument}}$：仪器精度误差（来自测量设备规格）
  \begin{itemize}
  \item CO$_2$ 连续分析仪：$\sim$0.1 ppm
  \item 烧瓶采样分析：$\sim$0.05 ppm
  \item 卫星 OCO-2：$\sim$1--2 ppm（柱浓度）
  \end{itemize}
\item $\sigma_{\text{representation}}$：代表性误差（观测点与模式网格的时空尺度差异）
  \begin{equation}
  \sigma_{\text{rep}} \approx \sqrt{\frac{\text{Var}(\mathbf{x}_{\text{subgrid}})}{\text{足迹权重}}},
  \end{equation}
  通常为仪器误差的 1--3 倍
\item $\sigma_{\text{transport}}$：输送模型误差（FLEXPART 不确定性）
  \begin{itemize}
  \item 边界层观测：10--20\% 模拟浓度
  \item 自由对流层观测：5--10\% 模拟浓度
  \end{itemize}
\end{itemize}

\subsubsection{自适应误差调节}

根据观测特征调整误差：
\begin{itemize}
\item \textbf{自由对流层观测}（$z > 3000\,\text{m}$）：$\sigma_i \leftarrow \beta_{\text{FT}} \cdot \sigma_i$（$\beta_{\text{FT}} > 1$，代表性误差增大）
\item \textbf{近海岸站点}（距海岸 $< 50\,\text{km}$）：$\sigma_i \leftarrow \beta_{\text{coast}} \cdot \sigma_i$（海陆通量混合增加不确定性）
\item \textbf{夜间观测}：$\sigma_i \leftarrow \beta_{\text{night}} \cdot \sigma_i$（边界层稳定，局地通量主导）
\end{itemize}

{\color{red}\textbf{[自适应调整的物理基础]：自由对流层空气经长时间混合，单点观测代表性强，但对地表通量敏感性降低，需放大误差以避免过拟合。海岸站点同时``看到''陆地和海洋通量，存在混合区不确定性。夜间边界层浅，输送模型对垂直混合参数化敏感，误差增大。}}

\subsection{输送算子 $\mathbf{H}$ 的构造}

输送算子 $\mathbf{H} \in \mathbb{R}^{M \times N}$ 连接地表通量与大气浓度，其物理意义为源-受体关系（Source-Receptor Relationship）：
\begin{equation}
H_{ij} = \frac{\partial y_i}{\partial x_j} = \text{第 $j$ 个网格单位通量对第 $i$ 个观测点混合比的贡献}.
\end{equation}

\subsubsection{FLEXPART 后向拉格朗日方法}

FLEXINVERT 通过 FLEXPART 后向轨迹计算 $\mathbf{H}$：
\begin{enumerate}
\item \textbf{粒子释放}：从观测点 $(x_{\text{obs}}, y_{\text{obs}}, z_{\text{obs}}, t_{\text{obs}})$ 释放 $J$ 个虚拟粒子（典型 $J \sim 10^4$）
\item \textbf{后向追踪}：使用风场数据（如 ERA5）向后积分粒子轨迹，通常追踪 10--30 天
\item \textbf{驻留时间统计}：记录粒子在每个网格 $(i,j,k)$ 和时间 $t$ 的驻留时间 $\Delta t'_{p,ijk,t}$
\item \textbf{密度归一化}：除以沿轨迹的空气密度 $\rho_{p,t}$
\item \textbf{垂直积分}：对边界层内的驻留时间积分
\item \textbf{粒子平均}：对所有粒子求平均
\end{enumerate}

数学表达式为：
\begin{equation}
H_{ij} = \frac{1}{J} \sum_{p=1}^{J} \sum_{k \in \text{PBL}} \frac{\Delta t'_{p,ij,k}}{\rho_{p,k}},
\end{equation}
单位为 (ppm$\cdot$m$^2\cdot$s)/kg，表示单位通量（kg/m$^2$/s）对混合比（ppm）的贡献。

{\color{red}\textbf{[后向轨迹的优势]：相比前向模拟（从所有源区释放粒子到观测点），后向模拟仅需从观测点释放，计算量随观测数 $M$ 而非网格数 $N$ 扩展。对于 $M \sim 10^3, N \sim 10^5$ 的典型问题，后向法效率提升 100 倍。这使高分辨率反演成为可能。}}

\subsubsection{足迹矩阵存储}

$\mathbf{H}$ 矩阵从不显式存储，而是作为``足迹文件''集合：
\begin{itemize}
\item 每个观测对应一个足迹文件：\texttt{footprint\_obs\_{i}.nc}
\item 文件包含：三维网格驻留时间场 $(n_x, n_y, n_t)$
\item 矩阵-向量乘积 $\mathbf{H}\mathbf{x}$ 通过卷积实现：
\begin{equation}
(\mathbf{H}\mathbf{x})_i = \sum_{j=1}^{N} H_{ij} x_j = \sum_{t,x,y} \text{footprint}_{i,txy} \cdot \text{flux}_{txy}.
\end{equation}
\item 伴随乘积 $\mathbf{H}^T\mathbf{y}$ 将观测空间梯度投影回通量空间
\end{itemize}

{\color{red}\textbf{[$\mathbf{H}$ 矩阵的稀疏性]：虽然 $\mathbf{H}$ 在数学上是 $M \times N$ 稠密矩阵，实际上极度稀疏。单个观测仅``看到''风向上游 1000--2000 km 范围的通量（$<1\%$ 网格），其余元素为零。这种隐式稀疏性通过足迹文件自然体现，无需显式存储零元素。}}

\subsubsection{域分解与背景贡献}

对于区域反演，输送算子分解为：
\begin{equation}
y_i^{\text{mod}} = \sum_{j \in \text{nest}} H_{ij}^{\text{nest}} x_j + \sum_{j \in \text{out}} H_{ij}^{\text{out}} x_j^{\text{fixed}} + H_i^{\text{bg}} y_i^{\text{bg}},
\end{equation}
其中：
\begin{itemize}
\item $\mathbf{H}^{\text{nest}}$：高分辨率嵌套域（待优化变量），如 $0.5^\circ \times 0.5^\circ$
\item $\mathbf{H}^{\text{out}}$：粗分辨率全球域（固定背景通量），如 $3^\circ \times 2^\circ$
\item $H^{\text{bg}}$：轨迹端点（10--30 天前）初始浓度贡献
\end{itemize}

背景浓度 $y_i^{\text{bg}}$ 通常来自全球模型（如 CarbonTracker, CAMS）。

\subsubsection{变分辨率网格聚合}

为降低计算成本，FLEXINVERT 采用投影算子 $\mathbf{\Gamma}$ 聚合低敏感度网格：
\begin{equation}
\mathbf{H}_{\text{vr}} = \mathbf{H}_{\text{hr}} \mathbf{\Gamma}^T,
\end{equation}
其中 $\mathbf{H}_{\text{hr}} \in \mathbb{R}^{M \times N_{\text{hr}}}$ 为高分辨率输送算子，$\mathbf{\Gamma} \in \mathbb{R}^{N_{\text{vr}} \times N_{\text{hr}}}$ 为聚合矩阵（$N_{\text{vr}} \ll N_{\text{hr}}$）。

聚合准则：
\begin{itemize}
\item 保留高足迹敏感度区域的细分辨率（$\sum_i H_{ij}^2 > \theta$）
\item 合并低敏感度邻近网格（$\sum_i H_{ij}^2 < \theta$）
\item 维持海陆边界分离（避免混合不同地表类型）
\end{itemize}

典型聚合可将 $N_{\text{hr}} \sim 10^6$ 降至 $N_{\text{vr}} \sim 10^4$--$10^5$，状态空间压缩 10--100 倍。

{\color{red}\textbf{[变分辨率的智慧]：观测网络通常集中在北美、欧洲等地区，这些``上游''区域需要细分辨率。而远离观测的区域（如南太平洋）对反演影响小，可粗化至 $5^\circ \times 5^\circ$ 甚至更大。这是``信息驱动的网格自适应''——分辨率随观测敏感度动态调整，既保证精度又控制计算量。}}

\subsection{数值求解策略}

FLEXINVERT 提供三种求解器，根据问题规模自动选择：

\subsubsection{解析解（小规模问题 $N < 10^4$）}

直接应用卡尔曼增益公式：
\begin{equation}
\mathbf{x}_a = \mathbf{x}_b + \mathbf{B}\mathbf{H}^T(\mathbf{H}\mathbf{B}\mathbf{H}^T + \mathbf{R})^{-1}(\mathbf{y}^{\text{obs}} - \mathbf{H}\mathbf{x}_b),
\end{equation}
需求逆 $M \times M$ 矩阵，复杂度 $\mathcal{O}(M^3 + M^2N)$。

后验协方差精确计算：
\begin{equation}
\mathbf{A} = (\mathbf{B}^{-1} + \mathbf{H}^T\mathbf{R}^{-1}\mathbf{H})^{-1} = \mathbf{B} - \mathbf{B}\mathbf{H}^T(\mathbf{H}\mathbf{B}\mathbf{H}^T + \mathbf{R})^{-1}\mathbf{H}\mathbf{B}.
\end{equation}

\subsubsection{共轭梯度法（中等规模 $10^4 < N < 10^6$）}

在 $\chi$ 空间迭代最小化 $J(\chi)$：
\begin{equation}
J(\chi) = \frac{1}{2}\chi^T\chi + \frac{1}{2}(\mathbf{H}\mathbf{B}^{1/2}\chi - \mathbf{d})^T\mathbf{R}^{-1}(\mathbf{H}\mathbf{B}^{1/2}\chi - \mathbf{d}),
\end{equation}
其中 $\mathbf{d} = \mathbf{y}^{\text{obs}} - \mathbf{H}\mathbf{x}_b$。

梯度计算：
\begin{equation}
\nabla_\chi J = \chi + \mathbf{B}^{1/2}\mathbf{H}^T\mathbf{R}^{-1}(\mathbf{H}\mathbf{B}^{1/2}\chi - \mathbf{d}).
\end{equation}

预条件共轭梯度（PCG）算法：
\begin{enumerate}
\item 初始化：$\chi^{(0)} = 0$，$\mathbf{r}^{(0)} = -\nabla J(\chi^{(0)})$，$\mathbf{p}^{(0)} = \mathbf{r}^{(0)}$
\item 迭代（$k = 0, 1, 2, \ldots$）：
\begin{align}
\alpha^{(k)} &= \frac{(\mathbf{r}^{(k)})^T\mathbf{r}^{(k)}}{(\mathbf{p}^{(k)})^T\nabla^2 J \mathbf{p}^{(k)}}, \\
\chi^{(k+1)} &= \chi^{(k)} + \alpha^{(k)} \mathbf{p}^{(k)}, \\
\mathbf{r}^{(k+1)} &= \mathbf{r}^{(k)} - \alpha^{(k)} \nabla^2 J \mathbf{p}^{(k)}, \\
\beta^{(k+1)} &= \frac{(\mathbf{r}^{(k+1)})^T\mathbf{r}^{(k+1)}}{(\mathbf{r}^{(k)})^T\mathbf{r}^{(k)}}, \\
\mathbf{p}^{(k+1)} &= \mathbf{r}^{(k+1)} + \beta^{(k+1)} \mathbf{p}^{(k)}.
\end{align}
\item 收敛准则：$\|\mathbf{r}^{(k)}\| < \epsilon \|\mathbf{r}^{(0)}\|$（$\epsilon \sim 10^{-6}$）
\end{enumerate}

\subsubsection{拟牛顿法 M1QN3（大规模 $N \sim 10^6$）}

有限内存 BFGS 方法，近似 Hessian 矩阵的逆：
\begin{equation}
\mathbf{H}_k \approx (\nabla^2 J)^{-1},
\end{equation}
通过最近 $m$ 次迭代的梯度差和步长更新（$m \sim 5$--20）。

M1QN3 算法特点：
\begin{itemize}
\item 仅需梯度计算，无需 Hessian
\item 内存需求 $\mathcal{O}(mN)$，远小于完整 BFGS 的 $\mathcal{O}(N^2)$
\item 超线性收敛速度（优于 PCG 的线性收敛）
\item 适用于光滑非线性问题（如 MCMC 初始化）
\end{itemize}

\subsubsection{求解器选择流程}

\begin{verbatim}
if N * M < 10^8:
    使用解析解（精确后验协方差）
elif N < 5 × 10^5 且问题条件数良好:
    使用 PCG（10--50 次迭代）
else:
    使用 M1QN3（对非线性更鲁棒）
\end{verbatim}

\subsubsection{计算复杂度对比}

\begin{table}[ht]
\centering
\caption{求解器复杂度对比（$N$ 为状态维度，$M$ 为观测数，$K$ 为迭代次数）}
\begin{tabular}{llll}
\hline
求解器 & 单次迭代 & 总复杂度 & 内存需求 \\
\hline
解析解 & - & $\mathcal{O}(M^3 + M^2N)$ & $\mathcal{O}(M^2 + MN)$ \\
PCG & $\mathcal{O}(MN)$ & $\mathcal{O}(K \cdot MN)$ & $\mathcal{O}(N + M)$ \\
M1QN3 & $\mathcal{O}(MN)$ & $\mathcal{O}(K \cdot MN)$ & $\mathcal{O}(mN)$ \\
\hline
\end{tabular}
\end{table}

典型问题规模特征：
\begin{itemize}
\item 状态向量维度：$N = 10^4$--$10^6$（空间网格 $\times$ 时间步）
\item 观测数量：$M = 10^2$--$10^4$（稀疏观测网络）
\item 条件数：$\kappa(\mathbf{B}) \sim 10^6$（需预条件化）
\item PCG 迭代次数：$K \sim 10$--100
\item M1QN3 迭代次数：$K \sim 20$--200
\item 计算时间：数分钟至数小时（取决于求解器与问题规模）
\end{itemize}

{\color{red}\textbf{[求解器权衡]：解析解虽然精确但仅适用于小问题（$N < 10^4$）。PCG 利用 $\chi$ 空间预条件化，将条件数从 $10^6$ 降至接近 1，收敛快。M1QN3 对初始点和非线性更鲁棒，常用于 MCMC 初始化或弱约束 4D-Var。实践中，先用 PCG 快速收敛至近似解，再用 M1QN3 精细调整，结合两者优势。}}

\end{document}
